#pragma once
#include <string>
#include "Grid.h"

using namespace std;

double omega = (double)2/3;

class MatrixVectorOperator
{
private:
public:
    MatrixVectorOperator(){};
    double max_norm(vector<double> v1, vector<double> v2)
    {
	if (v1.size() != v2.size())
	{
	    cout << "Error: unsupported size." << endl;
	    exit(0);
	}
	double maxn = 0;
	for (int i = 0; i < v1.size(); i++)
	{
	    if (v1[i] - v2[i] > maxn)
		maxn = v1[i] - v2[i];
	    else if (v2[i] - v1[i] > maxn)
		maxn = v2[i] - v1[i];
	}
	return maxn;
    };
    vector<vector<double>> get_A(int N)
    {
	vector<vector<double>> A(N-1);
	for (int i = 0; i < N-1; i++)
	{
	    A[i].resize(N-1);
	    for (int j = 0; j < N-1; j++)
	    {
		if (i == j)
		    A[i][j] = 2*(N*N);
		if (i == j+1 || j == i+1)
		    A[i][j] = -1*(N*N);
	    }
	}
	return A;
    };
    vector<vector<double>> get_D(int N)
    {
	vector<vector<double>> D(N-1);
        for (int i = 0; i < N-1; i++)
	{
	    D[i].resize(N-1);
	    D[i][i] = 0.5/(N*N);
	}
	return D;
    };
    vector<vector<double>> get_Tw(int N)
    {
	vector<vector<double>> Tw = get_A(N);
	for (int i = 0; i < N-1; i++)
	{
	    for (int j = 0; j < N-1; j++)
	    {
		if (i == j)
		    Tw[i][j] = 1-omega;
		if (i == j+1 || j == i+1)
		    Tw[i][j] = omega/2;
	    }
	}
	return Tw;
    };
    vector<vector<double>> get_Id(int N)
    {
	vector<vector<double>> Id(N/2 - 1);
	Id[0].resize(N - 1);
	Id[0][0] = 0.25;
	Id[0][1] = 0.5;
	Id[0][2] = 0.25;
	for (int i = 1; i < N/2 - 1; i ++)
	{
	    Id[i].resize(N - 1);
	    for (int j = 2; j < N - 1; j ++)
	    {
		Id[i][j] = Id[i-1][j-2];
	    }
	}
	return Id;
    };
    vector<vector<double>> get_Iu(int N)
    {
	vector<vector<double>> v = M_trans(get_Id(2*N));
	for (int i = 0; i < v.size(); i++)
	{
	    for (int j = 0; j < v[0].size(); j++)
	    {
		v[i][j] = 2.0*v[i][j];
	    }
	}
	return v;
    };
    vector<double> V_plus_V(vector<double> v1, vector<double> v2)
    {
	if (v1.size() != v2.size())
	{
	    cout << "Error: unsupported size." << endl;
	    exit(0);
	}
	vector<double> v(v1.size());
	for (int i = 0; i < v.size(); i++)
	{
	    v[i] = v1[i] + v2[i];
	}
	return v;
    };
    vector<double> V_minus_V(vector<double> v1, vector<double> v2)
    {
	if (v1.size() != v2.size())
	{
	    cout << "Error: unsupported size." << endl;
	    exit(0);
	}
	vector<double> v(v1.size());
	for (int i = 0; i < v.size(); i++)
	{
	    v[i] = v1[i] - v2[i];
	}
	return v;
    };
    vector<double> M_times_V(vector<vector<double>> m1, vector<double> v1)
    {
	if (v1.size() != m1[0].size())
	{
	    cout << "Matrix-vector product must conform!" << endl;
	    exit(0);
	}
	vector<double> v(m1.size());
	for (int i = 0; i < v.size(); i++)
	{
	    for (int j = 0; j < v1.size(); j++) {
                v[i] += m1[i][j] * v1[j];
            }
	}
	return v;
    };
    vector<vector<double>> C_times_M(double C, vector<vector<double>> m1)
    {
	int r = m1.size(), c = m1[0].size();
	vector<vector<double>> m(c);
	for (int i = 0; i < c; i++)
	{
	    m[i].resize(r);
	    for (int j = 0; j < r; j++)
	    {
		m[i][j] = m1[i][j]*C;
	    }
	}
	return m;
    };
    vector<double> C_times_V(double C, vector<double> v1)
    {
	int r = v1.size();
	vector<double> v(r);
	for (int i = 0; i < r; i++)
	{
	    v[i] = C*v1[i];
	}
	return v;
    };
    vector<vector<double>> M_trans(vector<vector<double>> m1)
    {
	int r = m1.size(), c = m1[0].size();
	vector<vector<double>> m(c);
	for (int i = 0; i < c; i++)
	{
	    m[i].resize(r);
	    for (int j = 0; j < r; j++)
	    {
		m[i][j] = m1[j][i];
	    }
	}
	return m;
    };
};
